Introduction à la modélisation statistique bayésienne

PRINCIPES DE L'ANALYSE DE DONNÉES BAYESIENNE

• Un ensemble de données / mesures observées à analyser :
  \( d_1, d_2, \dots, d_N \)

• Un modèle génératif défini par un ensemble de paramètres
  \( \theta_1, \theta_2, \dots, \theta_K \)

• Une connaissance a priori de ces paramètres

Nous avons les observations

The Bayesian approach places uncertainty not in the observations but rather in one’s lack of knowledge. For a Bayesian, the observed data are not uncertain—you observed what you observed.

Feinberg and Gonzalez (2012)

PRINCIPES DE L'ANALYSE DE DONNÉES BAYESIENNE

• Un ensemble de données / mesures observées à analyser :
  \( d_1, d_2, \dots, d_N \)

• Un modèle génératif défini par un ensemble de paramètres
  \( \theta_1, \theta_2, \dots, \theta_K \)

• Une connaissance a priori de ces paramètres

Nous cherchons

UTILISATION DE LA RÈGLE DE BAYES

Ces observations et ces paramètres sont regroupés dans un même modèle bayésien :

la distribution conjointe \( P(d_{1:N}~\theta_{1:K}) \).

2 possibilités pour décomposer cette distribution à l'aide du théorème de Bayes :

• En faisant apparaître la fonction de vraissemblance (likelihood) et la distribution a priori (prior)

\[ P(d_{1:N}~\theta_{1:K})~=~\color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})} \]

• En faisant apparaître la distribution postérieure

\[ P(d_{1:N}~\theta_{1:K})~=~\color{purple}{P(\theta_{1:K}~|~d_{1:N})}~\color{green}{P(d_{1:N})} \]

Ces deux décompositions sont égales !!!

INFÉRENCE BAYESIENNE

On infére les nouvelles valeurs des paramètres à partir des données

\[ \color{purple}{P(\theta_{1:K}~|~d_{1:N})} = \frac{\color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})}}{\color{green}{P(d_{1:N})}} \]

Ou encore

\[ \color{purple}{P(\theta_{1:K}~|~d_{1:N})} \propto \color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})} \]

L'objectif de l'analyse de données bayésienne

Passer d'une connaissance a priori sur les paramètres \( \color{steelblue}{P(\theta_{1:K})} \)
À une connaissance à postériori \( \color{purple}{P(\theta_{1:K}~|~d_{1:N})} \)

ANALYSE DE DONNÉES BAYÉSIENNE

1. Définir le modèle

• Identifier les paramètres du modèle génératif (Likelihood) associés aux données
  
• Définir une distribution a priori (prior) sur ces paramètres

2. Mettre à jour le modèle

• Calculer la distribution a posteriori des paramètres (ou une bonne approximation)

3. Interpréter la distribution postérieure (répondre à la question posée)

• Estimation du paramètre
• Comparaison de modèles

4. Vérification prédictive du modèle

• S'assurer que le modèle génère des données cohérentes avec les données observées

POURQUOI COMMENCER PAR CE MODÈLE

• Le modèle beta-binomial couvre un grand nombre de problèmes de la vie courante :

  • Réussite / Echec à un test psychométrique ?
  • Présence / absence d'effets secondaires lors d'un test d'un médicament ?
  • Résultat d'un questionnaire à réponse binaire vrai/faux
  • Estimation des résultats du deuxième tour de l'éléction présidentielle
  • …

• C'est un modèle simple

  • Un unique paramètre
  • Solution analytique

FIL ROUGE - lancer d'une pièce de monnaie

• Génération des données, la valeur de \( \theta \) étant connue (ici, \( \theta=0.6 \)).

library(tidyverse)
set.seed(1234567)
sample(x = c(0, 1), size = 50, replace = TRUE, prob = c(0.4,0.6)) %>%
        data.frame %>%
        mutate(x = seq_along(.), y = cumsum(.) / seq_along(.) ) %>%
        ggplot(aes(x = x, y = y), log = "y") +  geom_line(lwd = 1, col = "steelBlue") +
        geom_hline(yintercept = 0.5, lty = 3) +
        xlab("nombre de lancers") + ylab("proportion de faces") +
        ylim(0, 1) + theme_bw(base_size = 20)

LIKELIHOOD

• Nous considérons \( z \) comme étant le nombre de succès \( \{y_i\} \)
• Nous considérons N (le nombre d'observations) comme étant une constante
• Nous considérons \( \theta \) comme étant un paramètre de notre modèle

La fonction de likelihood s'écrit :

\[ L(\theta~|~z~N)=P(z~|~\theta~N)=\theta^z~(1-\theta)^{(N-z)} \]

c'est une fonction de la variable \( \theta \)

LIKELIHOOD - distribution Binomiale

Exemple :
Représentation graphique de la fonction de vraissemblance de \( \theta \) pour 10 lancers et 7 Faces

NbrTrial     = 100
N            = 10
nbrSuccess   = 7
theta        = seq(from = 0, to = 1, by = 1/NbrTrial)
likelihood   = (theta^(nbrSuccess)) * (1 - theta)^(N - nbrSuccess) 
data.frame(theta, likelihood) %>%
ggplot(aes(x = theta, y = likelihood)) + 
    geom_area(color="orangered", fill = "orangered", alpha =.4) +
    theme_bw(base_size = 20) + theme(legend.position="none",
          text = element_text(size=30))

PRIOR

Comment définir un prior dans le cas du lancer de pièce ?

Aspect sémantique, rendre compte ou modéliser :

• Une absence d'information
• La connaissance d'observations antérieures concernant la pièce lancée.
• Un niveau de certitude concernant ces observations antérieures

Aspect mathématique, pour une solution entièrement analytique :

• Les distribution postérieure et prior doivent avoir la même forme
• La vraissemblance marginale doit se calculer analytiquement

PRIOR - distribution Beta

Un candidat : La distribution beta

La distribution beta:

\[ \begin{align} p(\theta~|~a,b) &= beta(\theta~|~a,b) \\ &= \theta^{(a-1)}(1-\theta)^{(b-1)}/B(a,b) \\ &\propto \theta^{(a-1)}(1-\theta)^{(b-1)} \end{align} \]

où \( a \) et \( b \) sont deux paramètres tels que \( a\geq0~et~b\geq0 \) et \( B(a,b) \) est une constante de normalisation

Attention : Ne pas confondre dans R : la fonction de \( \theta \) et la constante (indépendante de \( \theta \))

a     = 2
b     = 3
theta = seq( from= 0 , to= 1 , length.out= 7 )
dbeta(theta,a,b)

[1] 0.0000000 1.3888889 1.7777778 1.5000000 0.8888889 0.2777778 0.0000000

beta(a,b)

[1] 0.08333333

PRIOR - distribution Beta

Vérifions les contraintes mathématiques

Le résultat du produit d'une distribution beta et d'une fonction de vraissemblance de Bernoulli est proportionnel à une distribution beta.
On dit que la distribution beta est la distribution conjugée de la fonction de vraissemblance de Bernoulli.

Soit un prior définit par \( ~~~\color{steelblue}{p(\theta~|~a,b) \propto \theta^{(a-1)}(1-\theta)^{(b-1)}} \)
Soit une fonction de vraissemblance associée à z “Face” pour N lancers \( ~~\color{orangered}{P(\{y_i\}~|~\theta)=\theta^z~(1-\theta)^{(N-z)}} \)

\[ \begin{align} \color{purple}{p(\theta~|~\{y_i\})} &\propto \color{orangered}{P(\{y_i\}~|~\theta)}~\color{steelblue}{p(\theta~|~a,b)} &&\mbox{inférence Bayésienne} \\ &\propto \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}} &&\mbox{application des formules précédentes} \\ &\propto \color{purple}{\theta^{(z+a-1)}~(1-\theta)^{(N-z+b-1)}} &&\mbox{en regroupant les termes identiques} \\ &\propto \color{purple}{\theta^{(A-1)}(1-\theta)^{(B-1)}} &&\mbox{avec } A=z+a \mbox{ et } B=N-z+b \end{align} \]

Conclusions : la distribution a priori et la distribution a posteriori ont la même forme

PRIOR - distribution Beta

Vérifions les aspects sémantiques

• On peut exprimer l'absence de connaissance à priori - a=1 et b=1 (courbe bleue).
• On modélise un biais en faveur de Face - a > b (courbe orange).
• On modélise un biais en faveur de Pile - a < b (courbe verte).

PRIOR - distribution Beta

Vérifions les aspects sémantiques

• Le niveau de certitude augmente avec la somme (\( \kappa = a+b \))

  • Aucune idée sur la provenance de la pièce : a = b = 1 -> prior faible (courbe bleue)
  • En attendant le début de l'expérience, on a lancé la pièce 10 fois pour 5 “Face” : a = b = 5 (courbe orange)
  • La pièce provient de la banque de France : a = b = 50 -> prior fort (courbe verte)

PRIOR - distribution Beta

Le plus souvent, on dispose :

• D'information sur la tendance centrale du paramètre (mode, moyenne)
• D'un niveau de certitude ( \( \kappa = a+b \) )

PRIOR - distribution Beta

Le mode
Supposons que l'on ait une estimation de la valeur la plus probable \( \omega \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior

\[ \begin{align} a &= \omega(\kappa-2)+1 \\ b &= (1-\omega)(\kappa-2)+1 &&\mbox{pour } \kappa>2 \end{align} \]

Pour \( ~~\omega=0.65 \) et \( ~~\kappa=25 \) alors \( ~~~p(\theta)=beta(\theta~|~15.95, 9.05)~ \) (courbe orange)
Pour \( ~~\omega=0.65 \) et \( ~~\kappa=10 \) alors \( ~~~p(\theta)=beta(\theta~|~6.2, 3.8)~ \) (courbe bleue)

PRIOR - distribution Beta

La moyenne
Supposons que l'on ait une estimation de la valeur moyenne \( \mu \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior

\[ \begin{align} a &= \mu \kappa \\ b &= (1-\mu) \kappa &&\mbox{ pour } \kappa>2 \end{align} \]

Pour \( \mu=0.8 \) et \( \kappa=7 \) alors \( p(\theta)=beta(\theta~|~5.6,1.4) \) (courbe orange)
Pour \( \mu=0.8 \) et \( \kappa=22 \) alors \( p(\theta)=beta(\theta~|~17.6,4.4) \) (courbe bleue)

PRIOR - distribution Beta

La moyenne et la variance
Supposons que l'on ait une estimation a priori de la valeur moyenne \( \mu \) et de la variance \( \sigma \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior

\[ \begin{align} a &= \mu~(\mu~(1-\mu)~/~\sigma^2-1) \\ b &= (1-\mu)~\mu~(\mu~(1-\mu)~/~\sigma^2-1) ,~~~~~ \mbox{ pour } \sigma < 0.28867 \\ \end{align} \]

Pour \( ~~\mu=0.8 \) et \( ~~\sigma=0.15 \) alors \( ~~~p(\theta)=beta(\theta~|~5.92, 3.19) \) (orange)
Pour \( ~~\mu=0.8 \) et \( ~~\sigma=0.1 \) alors \( ~~~p(\theta)=beta(\theta~|~14.14, 7.61) \) (bleue)

PRIOR - Son influence sur la moyenne de la distribution postérieure

Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :

\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]

Cas \( \mbox{ } N < a+b \)

PRIOR - Son influence sur la moyenne de la distribution postérieure

Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :

\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]

Cas \( \mbox{ } N = a+b \)

PRIOR - Son influence sur la moyenne de la distribution postérieure

Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :

\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]

Cas \( \mbox{ } N > a+b \)

CONCLUSION

There are two main desiderata for a mathematical description of data.

• First, the mathematical form should be comprehensible with meaningful parameters.
  
• The second desideratum for a mathematical description is that it should be descriptively adequate, which means, loosely, that the mathematical form should “look like” the data.

Kruschke (2015)

INTRODUCTION

The posterior distribution is always a compromise between the prior distribution and the likelihood function.
Kruschke (2015)

PRINCIPE

\[ Posterior = \frac{Likelihood~\times~Prior}{\color{orangered}{Marginal~Likelhood}} \]

PROBLÈME :

le calcul de la “marginal likelihood”, \( \color{green}{~p(D)} \)

\[ \begin{align} \color{green}{p(D)} &= \int P(D,\theta) \, \mathrm d\theta &&\mbox{marginalisation sur le paramètre} \theta \\ \color{green}{p(D)} &= \color{green}{\int P(D~|~\theta)~P(\theta) \, \mathrm d\theta} && \mbox{Application de la règle du produit} \end{align} \]

Dans la littérature, on le trouve sous le nom de : “evidence”, “probabilité des données” ou “Marginal Likelihood”

PRINCIPE

\[ Posterior = \frac{Likelihood~\times~Prior}{\color{orangered}{Marginal~Likelhood}} \]

Trois méthodes pour résoudre ce problème :

1. Solution analytique \( ~~\longrightarrow~~ \) Utilisation d'une likelihood conjuguée du prior
   
2. Solution discrètisée \( ~~\longrightarrow~~ \) Calcul de la solution sur un ensemble fini de points (Méthode Grid)
3. Solution approchée \( ~~\longrightarrow~~ \) On échantillonne “intelligemment” sur l'espace conjoint des paramètres (Méthode MCMC)

CAS ANALYTIQUE

Utilisation d'une likelihood conjuguée du prior : cas Beta-binomiale

\[ \begin{align} \color{purple}{p(\theta~|~z,N)} &= \color{orangered}{P(z~|~\theta,N)}~\color{steelblue}{p(\theta)/\color{green}{P(z~|~N)}} &&\mbox{Équation de Bayes} \\ &= \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\frac{\theta^{(a-1)}(1-\theta)^{(b-1)}}{B(a,b)}}/\color{green}{P(z~|~N)} &&\mbox{Application des formules précédetes} \\ &= \frac{\color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}}}{\color{steelblue}{B(a,b)}\color{green}{P(z~|~N)}} \\ &= \frac{\theta^{(\color{orangered}{z}+\color{steelblue}{a-1})}~(1-\theta)^{(\color{orangered}{N-z}+\color{steelblue}{b-1})}}{\color{steelblue}{B(a,b)}\color{green}{P(z~|~N)}} &&\mbox{En regroupant les termes identiques} \\ &= \frac{\theta^{(\color{orangered}{z}+\color{steelblue}{a-1})}~(1-\theta)^{(\color{orangered}{N-z}+\color{steelblue}{b-1})}}{B(\color{orangered}{z}+\color{steelblue}{a},\color{orangered}{N-z}+\color{steelblue}{b})} &&\mbox{Par identification à une fonction beta} \\ \end{align} \]

CAS ANALYTIQUE

Distributions discrètes

CAS ANALYTIQUE

Distributions continues

Le modèle (likelihood + prior) doit être défini à partir de l'interprétation que l'on peut faire des paramètres et des distributions
et non pour faciliter les calculs

MÉTHODE GRID

1. Définir la grille
   
2. Calculer la valeur du prior pour chaque valeur de la grille
   
3. Calculer la valeur de la likelihood pour chaque valeur de la grille
   
4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat

MÉTHODE GRID

1. Définir la grille
   
2. Calculer la valeur du prior pour chaque valeur de la grille
   
3. Calculer la valeur de la likelihood pour chaque valeur de la grille
   
4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat

MÉTHODE GRID

1. Définir la grille
   
2. Calculer la valeur du prior pour chaque valeur de la grille
   
3. Calculer la valeur de la likelihood pour chaque valeur de la grille
   
4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat

MÉTHODE GRID

1. Définir la grille
   
2. Calculer la valeur du prior pour chaque valeur de la grille
   
3. Calculer la valeur de la likelihood pour chaque valeur de la grille
   
4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat

MÉTHODE GRID - grille affinée

1. Définir la grille
   
2. Calculer la valeur du prior pour chaque valeur de la grille
   
3. Calculer la valeur de la likelihood pour chaque valeur de la grille
   
4. Calculer le produit prior x likelihood pour chaque valeur de la grille, puis normalisation du résultat

MÉTHODE GRID - Problème du nombre de paramètres

En affinant la grille on augmente le temps de calcul

• 3 paramètres avec une grille de \( 10^3 \) noeuds => une grille de \( 10^9 \) points de calcul
  
• 10 paramètres avec une grille de \( 10^3 \) noeuds => une grille de \( 10^{30} \) points de calcul
  

Le “superordinateur” chinois Tianhe-2 réalise \( 33,8 x 10^{15} \) opérations par seconde
Si on considère qu'il réalise 1 opération par noeud de la grille, il lui faudrait \( 2.96 x 10^{13} \) pour parcourir la grille une fois !!!
Par comparaison, l'âge de l'univers est approximativement de \( (4,354 ± 0,012)×10^{17} \) secondes

MÉTHODE PAR ÉCHANTILLONNAGE - Le principe

• On ne calcule pas explicitement la distribution postérieure
• On en donne une estimation à partir du produit likelihood x prior calculée sur un échantillon des noeuds de la grille

En pratique :

• La fonction sample de R

trajLength     = 250 # nbr of measure                    
theta          = 1:7
ptheta         = theta
trajectory     = sample(theta, prob = ptheta, size = trajLength, replace = TRUE)

• Les méthodes MCMC, il existe différentes implémentations.
  (Metropolis, Gibbs, Hamilton, …)

MÉTHODE PAR ÉCHANTILLONNAGE - precision

EN PRATIQUE

• Cas analytique : Une fonction

p_grid    =  seq( from=0 , to=1 , length.out=1000 )
a = b     = 1
N         = 9
z         = 6
posterior = dbeta(p_grid, z+a-1, N-z+b-1)

• Cas Discret : Une valeur par noeud de la grille

p_grid      = seq( from=0 , to=1 , length.out = 1000 )
prior       = rep( 1 , 1000 )
likelihood  = dbinom( 6 , size=8 , prob=p_grid )
posterior   = likelihood * prior
posterior   = posterior / sum(posterior)

• Cas échantillons : Un ensemble de valeurs par noeud de la grille

samples = sample( p_grid , prob = posterior , size = 1e4 , replace = TRUE )

CONCLUSION

Méthode analytique

• La distribution postérieure est décrite explicitement
• Le modèle est fortement contraint

Méthode Grid

• La distribution postérieure n'est donnée que pour un ensemble finie de valeurs
• Plus la grille est fine, meilleure est l'estimation de la distribution postérieure
• Compromis Précision - Temps de calcul => Faible nombre de paramètres

Méthodes par échantillonnage (MCMC)

• La distribution postérieure est décrite par un échantillon
• Vérification résultat

ESTIMATION DE LA TENDANCE CENTRALE

DESCRIPTION : A partir d'un échantillon d'une distribution postérieure, on peut calculer la moyenne, le mode et la médiane

library(rethinking)
chainmode( samples , adj=0.01 ) # pour le mode ou MAP
mean( samples )                 # pour la moyenne
median( samples )               # pour la médiane

Exemple pour un prior uniforme, et 10 lancers successifs et 3 Faces

ESTIMATION DE LA TENDANCE CENTRALE - Cas “pathologique”“

DESCRIPTION : A partir d'un échantillon d'une distribution postérieure, on peut calculer la moyenne, le mode et la médiane

library(rethinking)
chainmode( samples , adj=0.01 ) # pour le mode ou MAP
mean( samples )                 # pour la moyenne
median( samples )               # pour la médiane

Exemple pour un prior uniforme, et 3 lancers successifs et 3 Faces

ESTIMATION DE LA TENDANCE CENTRALE

DÉCISION : Comment chosir la valeur qui résume le “mieux” la distribution postérieure
Utiliser des fonctions d'espérance de gain ou, ce qui revient au même, de perte attendue

# pour la moyenne
sapply( p_grid , function(d) sum( posterior*( d - p_grid )^2 ) )
# pour la médiane
sapply( p_grid , function(d) sum( posterior*abs( d - p_grid ) ) )

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

• Le HPDI indique les valeurs du paramètre qui sont le plus crédibles
• Plus le HPDI est étroit et plus le degré de certitude est élevé
• La largeur du HPDI diminue avec l'augmentation du nombre de mesures
• Le pourcentage de la masse de probabilité n'est pas nécessairement 95

Définition: les valeurs du paramètre x contenues dans le HDI 89% sont telles que \( p(x)>W \) et

\[ \int_{x:p(x)>W} p(x) \, \mathrm dx=0.89 \]

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

En pratique - Package BEST

library(coda)
library(BEST)
set.seed(1234567)
iterMax = 10000
p_grid  = seq( from=0 , to=1 , length.out=iterMax )
pTheta  = dbeta( p_grid, 3, 10)
massVec = pTheta/sum(pTheta)
samples = sample( p_grid , size=1e4 , replace=TRUE , prob=pTheta )

#plotPost(samples, credMass = 0.85)
hdi( samples , credMass = 0.89 )

     lower      upper 
0.05280528 0.39613961 
attr(,"credMass")
[1] 0.89

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

En pratique - Package BEST

HIGHEST POSTERIOR DENSITY INTERVAL (HPDI)

En pratique - Package rethinking

library(rethinking)
set.seed(1234567)
iterMax = 10000
p_grid  = seq( from=0 , to=1 , length.out=iterMax )
pTheta  = dbeta( p_grid, 3, 10)
massVec = pTheta / sum(pTheta)
samples = sample( p_grid , size = 1e4 , replace = TRUE , prob = pTheta )

# utilise le package CODA pour calculer HPDI
HPDI( samples, prob = 0.89)

     |0.89      0.89| 
0.05280528 0.39613961 

ROPE - Region Of Practical Equivalence

On l'utilise pour “tester une hypothèse”

• Pour tester H0, la décision est basée sur le HPDI et sur le ROPE
• La valeur du paramètre (ici 0) est rejetée si son ROPE exclus le 95% HPDI
• La valeur du paramètre (ici 0) est accepté si le 95% HDI est entièrement inclus dans son ROPE
• Si le HDI et le ROPE se chevauchent on ne peut rien dire…

ROPE - Region Of Practical Equivalence

COMPARAISON DE MODÈLES - Principes

On considère une pièce. On la lance 200 fois et on obtient 115 Faces.

QUESTION : La pièce est-elle biaisée?

Nous construisons deux modèles et essayons de savoir lequel est le plus probable

\[ \begin{eqnarray*} \left\{ \begin{array}{ll} M_0~:~X \sim \mathcal{B}(N,0.5) & \quad \text{la pièce n'est pas biaisée}\\ M_1~:~X \sim \mathcal{B}(N,\theta) & \quad \text{la pièce est biaisée} \end{array}\right. \end{eqnarray*} \]

COMPARAISON DE MODÈLES - Bayes Factor

Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.

\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\color{orangered}{\frac{P(D~|~M_0)}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]

Dans cet exemple, le calcul de \( P(D~|~M_0) \) est facile

\( P(D~|~M_0) = C_{200}^{115}~(0.5)^{115}~(1-0.5)^{(200-115)} = 0.005955 \)

COMPARAISON DE MODÈLES - Bayes Factor

Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.

\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\frac{P(D~|~M_0)}{\color{orangered}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]

Pour calculer de \( P(D~|~M_1) \), c'est plus compliqué : \( \color{orangered}{~~~\textbf{On ne connait pas }\theta} \)

Le calcul de la vraissemblance dans ce cas se fait par marginalisation sur le paramètre \( \theta \) \[ \begin{align} P(D~|~M_1) &= \int_0^1 P(D~|~M_1, \theta) P(\theta) \, \mathrm d\theta \\ &= \int_0^1 C^z_N~\theta^z~(1-\theta)^{N-z}~\text{d}\theta \\ &= N+1 \end{align} \]

COMPARAISON DE MODÈLES - Bayes Factor

Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.

\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\color{orangered}{\frac{P(D~|~M_0)}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]

Dans l'exemple,

\[ B_{01} = \frac{P(M_0~|~D)}{P(M_1~|~D)} = \frac{0.005955}{0.005} = 1.1971 \]

Le rapport de probabilités a augmenté de 20% après avoir pris connaissance des données.
La balance penche ici en faveur du modèle nul, \( B_{01} > 1 \)

COMPARAISON DE MODÈLES - Autre exemple

BIC critère d’information bayésien

Schwarz (1977) propose une approximation du log de la vraissemblance:

La vraissemblance du modèle pas toujours calculable : \[ P(D~|~M)=\int_{\theta} P(D~|~M, \theta) \, \mathrm d\theta \]

Calcul d'une valeur approchée \[ ln(P(D~|~\theta)) \approx ln(L(\widehat{\theta})) - \frac{t}{2}~ln(N) \]

avec \( L(\widehat{\theta}) \) le maximum de vraissemblance
\( t \) le nombre de paramètres
\( N \) le nombre d'observations

MODEL CHECKING

Les deux rôles de la fonction de vraissemblance

• C'est une fonction de \( \theta \) pour le calcul de la distribution postérieure : \( L(\theta~|~z~N) \)
• Lorsque \( \theta \) est fixé, c'est une distribution de probablilité : \( P(z~|~\theta~N)=\theta^z~(1-\theta)^{(N-z)} \)

Exemple
Générer 10000 valeurs à partir d'une loi binomiale basée sur 9 lancers et une probabilité de Face de 0,6 :

w <- rbinom(n = 10e4, size = 9, prob = 0.6)

MODEL CHECKING

Deux sources d'incertitude dans ces prédictions

• Incertitude lié à la valeur de \( p \) (évidemment !)
  -> Chaque valeur apparaît avec une probabilité \( p \)
  
• Incertitude sur la valeur de \( p \) elle-même
  -> Pour chaque valeur de \( p \) on peut calculer une distribution implicite

Exemple
Générer 10000 valeurs à partir d'une loi binomiale basé sur 9 lancers et une probabilité de Face de 0,6
pour chaque valeur de p

w <- rbinom(n = 10e4, size = 9, prob = samples)

MODEL CHECKING